#### libraries
library(xtable)
library(foreign)
library(dplyr)
library(ggplot2)
library(tidyverse)
library(reshape2)
library(texreg)

#### functions
mean.na <- function(x){
  mean(x,na.rm=T)
}

#### load data 
load("~/Dropbox/NAES/stable-replication-data-05062022.Rdata")

#### descriptive statistics table ----
dta22$VOTE16 <- (dta22$NONVOTE16==0)*1
dtasub <- dta22[! dta22$TRUMPVBIDEN15B %in% c(NA) & ! dta22$TRUMPVC11B %in% c(NA),]

vars <- c("EDYEARS6","FEMALE6","BLACK6","HISP6","WHITE6","INCOME6","UNION0","AGE6","PROTESTANT0","CATHOLIC0","VOTE16")

rmat <- matrix(NA,length(vars),6)
for(i in 1:length(vars)){
  txt <- paste("hold1 <- dta22$",vars[i],sep="")
  eval(parse(text=txt))
  rmat[i,1] <- mean(hold1,na.rm=T)
  rmat[i,2] <- sd(hold1,na.rm=T)
  rmat[i,3] <- length(na.omit(hold1))
  
  txt <- paste("hold <- dtasub$",vars[i],sep="")
  eval(parse(text=txt))
  rmat[i,4] <- mean(hold,na.rm=T)
  rmat[i,5] <- sd(hold,na.rm=T)
  rmat[i,6] <- length(na.omit(hold))
}

rownames(rmat) <- vars

xtable(rmat[,c(1,2,4,5)],digits=c(0,3,3,3,3))

#### ELECTION FAIRNESS FIGURE ----
idx <- which((! dta22$ELECTFAIR1 %in% c(NA)) & (! dta22$ELECTFAIR12  %in% c(NA))  & 
               (! dta22$ELECTFAIR15  %in% c(NA))) 
length(idx)

dta23 <- dta22[idx,]
dim(dta23)

rmat <- matrix(NA,3,5)
rmat[1,1] <- mean(dta23$ELECTFAIR1,na.rm=T)
rmat[1,2] <- mean(dta23$ELECTFAIR8,na.rm=T)
rmat[1,3] <- mean(dta23$ELECTFAIR9,na.rm=T)
rmat[1,4] <- mean(dta23$ELECTFAIR12,na.rm=T)
rmat[1,5] <-mean(dta23$ELECTFAIR15,na.rm=T)

rmat[2,1] <-mean(dta23$ELECTFAIR1[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,2] <-mean(dta23$ELECTFAIR8[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,3] <-mean(dta23$ELECTFAIR9[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,4] <- mean(dta23$ELECTFAIR12[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,5] <- mean(dta23$ELECTFAIR15[dta23$PID1 %in% c(1,2,3)],na.rm=T)

rmat[3,1] <- mean(dta23$ELECTFAIR1[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,2] <- mean(dta23$ELECTFAIR8[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,3] <- mean(dta23$ELECTFAIR9[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,4] <- mean(dta23$ELECTFAIR12[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,5] <- mean(dta23$ELECTFAIR15[dta23$PID1 %in% c(5,6,7)],na.rm=T)

#pdf("election-fairness-2008-03022021.pdf")
#jpeg("election-fairness-01082021.jpg")
plot(rmat[1,],x=c(2008,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="Year",pch=16,lty=1,type="b",ylab="Election Fairness",
     main="Perceived Election Fairness (1-5)",cex.main=1.8,cex.lab=1.3)
par(new=T)
plot(rmat[2,],x=c(2008,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="",pch=18,lty=1,type="b",col="blue",ylab="",xaxt="n",yaxt="n")
par(new=T)
plot(rmat[3,],x=c(2008,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="",pch=17,lty=1,type="b",col="red",ylab="",xaxt="n",yaxt="n")
text(x=2009.75,y=4.05,col="red","'08 Republicans")
text(x=2009.75,y=3.4,col="blue","'08 Democrats")
#dev.off()

#### rather live 
idx <- which((! dta22$LEGIT6A %in% c(NA)) & (! dta22$LEGIT11A  %in% c(NA))  & 
               (! dta22$LEGIT15A  %in% c(NA))) 
length(idx)

dta23 <- dta22[idx,]

rmat <- matrix(NA,3,5)
rmat[1,1] <- mean(dta23$LEGIT6A,na.rm=T)
rmat[1,2] <- mean(dta23$LEGIT8A,na.rm=T)
rmat[1,3] <- mean(dta23$LEGIT9A,na.rm=T)
rmat[1,4] <- mean(dta23$LEGIT11A,na.rm=T)
rmat[1,5] <-mean(dta23$LEGIT15A,na.rm=T)

rmat[2,1] <-mean(dta23$LEGIT6A[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,2] <-mean(dta23$LEGIT8A[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,3] <-mean(dta23$LEGIT9A[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,4] <-mean(dta23$LEGIT11A[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,5] <- mean(dta23$LEGIT15A[dta23$PID1 %in% c(1,2,3)],na.rm=T)

rmat[3,1] <- mean(dta23$LEGIT6A[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,2] <- mean(dta23$LEGIT8A[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,3] <- mean(dta23$LEGIT9A[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,4] <- mean(dta23$LEGIT11A[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,5] <- mean(dta23$LEGIT15A[dta23$PID1 %in% c(5,6,7)],na.rm=T)

### reverse polarity
rmat <- 6-rmat

#pdf("/users/danielhopkins/Dropbox/NAES/figures/legit-rather-live-2008-pid-03022021.pdf")
#jpeg("/users/danhop/Dropbox/NAES/legit-rather-live-2008-pid-01082021.jpg")
plot(rmat[1,],x=c(2012.83333,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="Year",pch=16,lty=1,type="b",ylab="Agree/Disagree",
     main="Rather Live in Our System",cex.main=1.8,cex.lab=1.3)
par(new=T)
plot(rmat[2,],x=c(2012.83333,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="",pch=18,lty=1,type="b",col="blue",ylab="",xaxt="n",yaxt="n")
par(new=T)
plot(rmat[3,],x=c(2012.83333,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="",pch=17,lty=1,type="b",col="red",ylab="",xaxt="n",yaxt="n")
text(x=2013.75,y=4.55,col="red","'08 Republicans")
text(x=2013.75,y=4.05,col="blue","'08 Democrats")
#dev.off()

### not critical
idx <- which((! dta22$LEGIT6D %in% c(NA)) & (! dta22$LEGIT11D  %in% c(NA))  & 
               (! dta22$LEGIT15D  %in% c(NA))) 
length(idx)
dta23 <- dta22[idx,]

rmat <- matrix(NA,3,5)
rmat[1,1] <- mean(dta23$LEGIT6D,na.rm=T)
rmat[1,2] <- mean(dta23$LEGIT8D,na.rm=T)
rmat[1,3] <- mean(dta23$LEGIT9D,na.rm=T)
rmat[1,4] <- mean(dta23$LEGIT11D,na.rm=T)
rmat[1,5] <-mean(dta23$LEGIT15D,na.rm=T)

rmat[2,1] <-mean(dta23$LEGIT6D[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,2] <-mean(dta23$LEGIT8D[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,3] <-mean(dta23$LEGIT9D[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,4] <-mean(dta23$LEGIT11D[dta23$PID1 %in% c(1,2,3)],na.rm=T)
rmat[2,5] <- mean(dta23$LEGIT15D[dta23$PID1 %in% c(1,2,3)],na.rm=T)

rmat[3,1] <- mean(dta23$LEGIT6D[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,2] <- mean(dta23$LEGIT8D[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,3] <- mean(dta23$LEGIT9D[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,4] <- mean(dta23$LEGIT11D[dta23$PID1 %in% c(5,6,7)],na.rm=T)
rmat[3,5] <- mean(dta23$LEGIT15D[dta23$PID1 %in% c(5,6,7)],na.rm=T)

# reverse polarity
rmat <- 6-rmat

#pdf("/users/danielhopkins/Dropbox/NAES/figures/legit-critical-2008-pid-03022021.pdf")
plot(rmat[1,],x=c(2012.83333,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="Year",pch=16,lty=1,type="b",ylab="Agree/Disagree",
     main="Not Critical of Our System",cex.main=1.8,cex.lab=1.3)
par(new=T)
plot(rmat[2,],x=c(2012.833333,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="",pch=18,lty=1,type="b",col="blue",ylab="",xaxt="n",yaxt="n")
par(new=T)
plot(rmat[3,],x=c(2012.833333,2014.8333,2014.916,2016.916,2020.83333),ylim=c(1,5),
     xlab="",pch=17,lty=1,type="b",col="red",ylab="",xaxt="n",yaxt="n")
text(x=2019.15,y=2.65,col="red","'08 Republicans")
text(x=2019.15,y=1.65,col="blue","'08 Democrats")
#dev.off()

##### feeling thermometer figures ----

dta22$TRUMPFT11 <- as.numeric(as.character(dta22$TRUMPFT11))

dta22.sub <- dta22[! dta22$TRUMPFT10 %in% c(NA) & ! dta22$TRUMPFT15 %in% c(NA),]

fttrump <- ggplot(dta22.sub, aes(x=TRUMPFT10)) + 
  geom_density(color="darkblue", fill="lightblue",alpha=.4)+
  geom_density(aes(x=TRUMPFT15),color="darkred",fill="red",alpha=.4)+
  labs(title="Feeling Thermometers: Trump",x="Rating", y = "Share")+
  geom_text(x=50, y=0.0085, label="Jan. 2016")+
  geom_text(x=90, y=0.011, label="Oct. 2020")

#pdf("~/Dropbox/NAES/figures/trumpft-02152021.pdf")
fttrump
#dev.off()

dta27.sub <- dta22[! dta22$TRUMPFT15 %in% c(NA) & ! dta22$FTOBAMA6 %in% c(NA),]
quantile(dta27.sub$TRUMPFT15,.9,na.rm=T)

fttrumpobama.re <- ggplot(dta27.sub, aes(x=FTOBAMA6)) + 
  geom_density(color="darkblue", fill="lightblue",alpha=.4)+
  geom_density(aes(x=TRUMPFT15),color="darkred",fill="red",alpha=.4)+
  labs(title="Feeling Thermometers at Re-Election: \n Trump and Obama",x="Rating", y = "Share")+
  geom_text(x=50, y=0.0045, label="Trump 2020")+
  geom_text(x=70, y=0.011, label="Obama 2012")

#pdf("~/Dropbox/NAES/figures/trump-obama-reelect.pdf")
fttrumpobama.re
#dev.off()

#### appendix feeling thermometer figures ----
dta23.sub <- dta22[! dta22$CLINTONFT11N %in% c(NA) & ! dta22$CLINTONFT3 %in% c(NA),]
dim(dta23.sub)
quantile(dta23.sub$CLINTONFT11N,.9)

ftclinton <- ggplot(dta23.sub, aes(x=CLINTONFT3)) + 
  geom_density(color="darkblue", fill="lightblue",alpha=.4)+
  geom_density(aes(x=CLINTONFT11N),color="darkred",fill="red",alpha=.4)+
  labs(title="Feeling Thermometers: Clinton",x="Rating", y = "Share")+
  geom_text(x=50, y=0.0135, label="Spring/Summer 2008")+
  geom_text(x=20, y=0.0115, label="Oct. 2016")

#pdf("clintonft-02152021.pdf")
ftclinton
#dev.off()

dta24.sub <- dta22[! dta22$OBAMAFT3 %in% c(NA) & ! dta22$FTOBAMA6 %in% c(NA),]
dim(dta24.sub)
quantile(dta24.sub$FTOBAMA6,.9)
quantile(dta24.sub$OBAMAFT13,.9,na.rm=T)

ftobama <- ggplot(dta24.sub, aes(x=OBAMAFT3)) + 
  geom_density(color="darkblue", fill="lightblue",alpha=.4)+
  geom_density(aes(x=FTOBAMA6),color="darkred",fill="red",alpha=.4)+
  labs(title="Feeling Thermometers: Obama",x="Rating", y = "Share")+
  geom_text(x=50, y=0.0125, label="Spring/Summer 2008")+
  geom_text(x=75, y=0.0135, label="Oct. 2012")

#pdf("obamaft-02152021.pdf")
ftobama
#dev.off()

dta25.sub <- dta22[! dta22$ROMNEYFT1 %in% c(NA) & ! dta22$FTROMNEY6 %in% c(NA),]
dim(dta25.sub)
quantile(dta22$FTROMNEY6,na.rm=T,.9)

ftromney <- ggplot(dta25.sub, aes(x=ROMNEYFT1)) + 
  geom_density(color="darkblue", fill="lightblue",alpha=.4)+
  geom_density(aes(x=FTROMNEY6),color="darkred",fill="red",alpha=.4)+
  labs(title="Feeling Thermometers: Romney",x="Rating", y = "Share")+
  geom_text(x=40, y=0.015, label="Winter 2008")+
  geom_text(x=82.5, y=0.012, label="Oct. 2012")

#pdf("romneyft-02152021.pdf")
ftromney
#dev.off()

#### 

sum(dta27.sub$TRUMPFT14 > 89,na.rm=T)/sum(dta27.sub$TRUMPFT14 %in% 0:100,na.rm=T)
sum(dta27.sub$FTOBAMA6 > 89,na.rm=T)/sum(dta27.sub$FTOBAMA6 %in% 0:100,na.rm=T)

sum(dta27.sub$TRUMPFT14 > 94,na.rm=T)/sum(dta27.sub$TRUMPFT14 %in% 0:100,na.rm=T)
sum(dta27.sub$FTOBAMA6 > 94,na.rm=T)/sum(dta27.sub$FTOBAMA6 %in% 0:100,na.rm=T)

M <- 1000
obama <- trump <- c()
for(i in 1:M){
  rs <- sample(1:dim(dta27.sub)[1],dim(dta27.sub)[1],replace=T)
  dta27.rs <- dta27.sub[rs,]
  trump[i] <- sum(dta27.rs$TRUMPFT15 > 94,na.rm=T)/sum(dta27.rs$TRUMPFT15 %in% 0:100,na.rm=T)
  obama[i] <- sum(dta27.rs$FTOBAMA6 > 94,na.rm=T)/sum(dta27.rs$FTOBAMA6 %in% 0:100,na.rm=T)
}

1-sum(trump > obama)/M

mean(obama)
sum(trump > obama)/length(trump < obama)


mean(dta27.sub$FTOBAMA6,na.rm=T)
mean(dta27.sub$TRUMPFT15,na.rm=T)

sd(dta27.sub$FTOBAMA6,na.rm=T)
sd(dta27.sub$TRUMPFT15,na.rm=T)

##### partisanship ----
dta23.sub <- dta22[! dta22$TRUMPVBIDEN15B %in% c(NA) & ! dta22$PID11 %in% c(NA) & ! dta22$PID4 %in% c(NA),]

tt <- table(dta23.sub$TRUMPVBIDEN15B,dta23.sub$PID11)
t2020 <- tt[3,]/apply(tt,2,sum)

tt2 <- table(dta23.sub$ROMNEYVO6B,dta23.sub$PID4)
t2012 <- tt2[3,]/apply(tt2,2,sum)

#pdf("~/Dropbox/NAES/figures/pid-lagged-party-02192021.pdf",width=10)
plot(x=seq(from=1,to=7),y=t2012,ylim=c(0,1),pch=16,type="b",
     xaxt="n",xlab="",ylab="Share Backing GOP",
     main="GOP Share by Lagged PID",
     cex.main=1.5,cex.lab=1.25)
par(new=T)
plot(x=seq(from=1,to=7),y=t2020,ylim=c(0,1),pch=17,type="b",
     xaxt="n",xlab="",ylab="",lty=2)
axis(side=1,at=seq(from=1,to=7),labels=c("Str D","2","3","Ind","5","6","Str R"))
text(x=4.5,y=.1,"2012")
text(x=4,y=.4,"2020")
dev.off()

#### partisan stability ----

#So, on the right-most end, it would be: 0.2930818 (3 x GOP)
#then, next on the right: 0.06792453 (2 x GOP, 1 Undecided)
#then: 0.02138365 (2 x Undecided, 1 x GOP)
#third from the left: 0.03522013 (1 x Dem, 2 Undecided)
#second from the left: 0.0754717 (2 x Dem, 1 Undecided)
#On the left-most end, it would be: 0.3647799 (3 x Dem)
vote_choices <- data.frame(vote_stability = c("3 x Dem", "2 x Dem, 1 Undecided", "1 x Dem, 2 Undecided", "", "1 x GOP, 2 Undecided", "2 x GOP, 1 Undecided", "3 x GOP"), 
                           Proportion = c(36, 8, 4, 14, 2, 7, 29))
vote_choices <- vote_choices %>% 
  mutate(Year = "2006")
vote_choices.melt <- melt(vote_choices)
vote_choices.melt$quality <- factor(vote_choices.melt$vote_stability, levels = vote_choices$vote_stability)
ggplot(vote_choices.melt, aes(x = variable, y = value, fill = fct_rev(quality))) +
  geom_col() +
  geom_text(aes(label = paste0(value, "%")),
            position = position_stack(vjust = 0.5)) +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal(base_size = 18) +
  ylab("Percentage") +
  xlab(NULL) + guides(fill=guide_legend(title="'12, '16, '20 votes")) +
  scale_fill_brewer(palette = "RdBu")


###### turnout & partisanship ----
#load("/users/danielhopkins/Dropbox/NAES/2020wave/merged-2020-n2471-wcatalist-12082021.Rdata")

#dta.all4$NONVOTE08 <- 1*(dta.all4$E2008GVM.x=="")
#dta.all4$NONVOTE10 <- 1*(dta.all4$E2010GVM.x=="")
#dta.all4$NONVOTE12 <- 1*(dta.all4$E2012GVM.x=="")
#dta.all4$NONVOTE14 <- 1*(dta.all4$E2014GVM.x=="")
#dta.all4$NONVOTE16 <- 1*(dta.all4$E2016GVM.x=="")

dta.dem <- dta22[dta22$PID6 %in% c(1,2,3),]
dta.gop <- dta22[dta22$PID6 %in% c(5,6,7),]
dta.ind <- dta22[dta22$PID6 %in% c(4),]
dta.all <- dta22[dta22$PID6 %in% c(1:7),]

gop12 <- mean(dta.gop$NONVOTE12)
gop14 <- mean(dta.gop$NONVOTE14)
gop16 <- mean(dta.gop$NONVOTE16)
gop18 <- mean(dta.gop$NONVOTE18,na.rm=T)

dem12 <- mean(dta.dem$NONVOTE12)
dem14 <- mean(dta.dem$NONVOTE14)
dem16 <- mean(dta.dem$NONVOTE16)
dem18 <- mean(dta.dem$NONVOTE18,na.rm=T)

all12 <- mean(dta.all$NONVOTE12)
all14 <- mean(dta.all$NONVOTE14)
all16 <- mean(dta.all$NONVOTE16)
all18 <- mean(dta.all$NONVOTE18,na.rm=T)

#pdf("~/Dropbox/NAES/2016wave/turnout-by-partisanship-12192021.pdf")

plot(c(1-gop12,1-gop14,1-gop16,1-gop18),x=c(2012,2014,2016,2018),type="b",pch=16,ylim=c(0,1),lwd=2,ylab="Share Voting",xlab="",xaxt="n",
     main="Validated Turnout by \n Partisanship: Panel Data",cex.main=1.8,cex.lab=1.3,xlim=c(2012,2020))
axis(side=1,at=c(2012,2014,2016,2018),labels=c(2012,2014,2016,2018))
par(new=T)
plot(c(1-dem12,1-dem14,1-dem16,1-dem18),x=c(2012,2014,2016,2018),type="b",pch=17,ylim=c(0,1),lwd=2,ylab="",xlab="",lty=2,xaxt="n",xlim=c(2012,2020))
par(new=T)
plot(c(1-all12,1-all14,1-all16,1-all18),x=c(2012,2014,2016,2018),type="b",pch=18,ylim=c(0,1),lwd=2,ylab="",xlab="",lty=3,xaxt="n",xlim=c(2012,2020))

text(x=2014,y=.79,"Republicans",cex=1.8)
text(x=2014,y=.55,"Democrats",cex=1.8)
text(x=2019,y=1-all18,"All",cex=1.8)

#dev.off()

#### ideological extremity

dta22.sub <- dta22[! dta22$IDEOTRUMP11 %in% c(NA) & 
                     ! dta22$IDEOTRUMP15 %in% c(NA),]

tt11 <- table(dta22.sub$IDEOTRUMP11-dta22.sub$IDEO10)/sum(table(dta22.sub$IDEOTRUMP11-dta22.sub$IDEO10))
tt15 <- table(dta22.sub$IDEOTRUMP15-dta22.sub$IDEO10)/sum(table(dta22.sub$IDEOTRUMP15-dta22.sub$IDEO10))

df11 <- rbind(tt11,tt15)

#pdf("~/Dropbox/NAES/figures/trump-relative-ideo-06182021.pdf",width=12)
bp <- barplot(df11,
              col=c("black","lightgray"),
              beside=T,
              main="Perceived Relative Trump Ideology",
              xlab="Trump More Liberal than Resp.                          Trump More Conservative than Resp.")
text(x=10.5,y=.07,"Oct. '16")
text(x=29.5,y=.15,"Oct. '20")
legend(x=3,y=.15,fill=c("black","lightgray"),
       c("Oct. '16; Mean=-0.72",
         "Oct. '20; Mean=-1.42"))
#dev.off()

tt11 <- table(dta22$IDEOBIDEN15-dta22$IDEO10)/sum(table(dta22$IDEOBIDEN15-dta22$IDEO10))
tt15 <- table(dta22$IDEOTRUMP15-dta22$IDEO10)/sum(table(dta22$IDEOTRUMP15-dta22$IDEO10))
tt11a <- c(tt11,0)

df12 <- rbind(tt11a,tt15)

#pdf("~/Dropbox/NAES/figures/trump-biden-ideo-06182021.pdf",width=12)
bp <- barplot(df12,
              col=c("black","lightgray"),
              beside=T,
              main="Perceived Relative Trump/Biden Ideology",
              xlab="Resp. More Liberal                               Resp. More Conservative")
text(x=10.5,y=.1325,"Trump")
text(x=29.5,y=.15,"Biden")
legend(x=3,y=.185,fill=c("black","lightgray"),c("Trump; Mean=-1.42",
                                                "Biden; Mean=2.10"))
dev.off()

tt11 <- table(dta22.sub$OBAMAIDEO07-dta22.sub$IDEO7)/sum(table(dta22.sub$OBAMAIDEO07-dta22.sub$IDEO7))
tt15 <- table(dta22.sub$ROMNEYIDEO07-dta22.sub$IDEO7)/sum(table(dta22.sub$ROMNEYIDEO07-dta22.sub$IDEO7))
#table(dta22.sub$IDEOCLINTON11-dta22.sub$IDEO10)/sum(table(dta22.sub$IDEOCLINTON11-dta22.sub$IDEO10))
df11 <- rbind(tt11,tt15)

mean(dta22.sub$OBAMAIDEO07-dta22.sub$IDEO7,na.rm=T)
mean(dta22.sub$ROMNEYIDEO07-dta22.sub$IDEO7,na.rm=T)

#pdf("~/Dropbox/NAES/figures/romney-obama-ideo-06182021.pdf",width=12)
bp <- barplot(df11,
              col=c("black","lightgray"),
              beside=T,
              main="Perceived Relative Romney/Obama Ideology",
              xlab="Resp. More Liberal                               Resp. More Conservative")
text(x=10.5,y=.1325,"Romney")
text(x=29.5,y=.15,"Obama")
legend(x=3,y=.185,fill=c("black","lightgray"),c("Romney; Mean=-1.38",
                                                "Obama; Mean=2.08"))
dev.off()


##### prejudice measures -----

dta.sub <- dta22[! dta22$WBPREJUDICE15 %in% c(NA),]

flag <- ! dta22$WBPREJUDICE11 %in% c(NA) & ! dta22$WBPREJUDICE13 %in% c(NA) & dta22$WHITE6==1
sum(flag)
d.flag <- ! dta22$WBPREJUDICE11 %in% c(NA) & ! dta22$WBPREJUDICE13 %in% c(NA) & dta22$PID4 %in% c(1,2,3)  & dta22$WHITE6==1
sum(d.flag)
r.flag <- ! dta22$WBPREJUDICE11 %in% c(NA) & ! dta22$WBPREJUDICE13 %in% c(NA) & dta22$PID4 %in% c(5,6,7)  & dta22$WHITE6==1
sum(r.flag)

prejudice.scale.all <- c(
  mean(dta22$WBPREJUDICE4[flag],na.rm=T),
  mean(dta22$WBPREJUDICE6[flag],na.rm=T),
  mean(dta22$WBPREJUDICE7[flag],na.rm=T),
  mean(dta22$WBPREJUDICE10[flag],na.rm=T),
  mean(dta22$WBPREJUDICE11[flag],na.rm=T),
  mean(dta22$WBPREJUDICE12[flag],na.rm=T),
  mean(dta22$WBPREJUDICE13[flag],na.rm=T),
  mean(dta22$WBPREJUDICE14[flag],na.rm=T),
  mean(dta22$WBPREJUDICE15[flag],na.rm=T)
)

prejudice.scale.d <- c(
  mean(dta22$WBPREJUDICE4[d.flag],na.rm=T),
  mean(dta22$WBPREJUDICE6[d.flag],na.rm=T),
  mean(dta22$WBPREJUDICE7[d.flag],na.rm=T),
  mean(dta22$WBPREJUDICE10[d.flag],na.rm=T),
  mean(dta22$WBPREJUDICE11[d.flag],na.rm=T),
  mean(dta22$WBPREJUDICE12[d.flag],na.rm=T),
  mean(dta22$WBPREJUDICE13[d.flag],na.rm=T),
  mean(dta22$WBPREJUDICE14[d.flag],na.rm=T),
  mean(dta22$WBPREJUDICE15[d.flag],na.rm=T)
)

prejudice.scale.r <- c(
  mean(dta22$WBPREJUDICE4[r.flag],na.rm=T),
  mean(dta22$WBPREJUDICE6[r.flag],na.rm=T),
  mean(dta22$WBPREJUDICE7[r.flag],na.rm=T),
  mean(dta22$WBPREJUDICE10[r.flag],na.rm=T),
  mean(dta22$WBPREJUDICE11[r.flag],na.rm=T),
  mean(dta22$WBPREJUDICE12[r.flag],na.rm=T),
  mean(dta22$WBPREJUDICE13[r.flag],na.rm=T),
  mean(dta22$WBPREJUDICE14[r.flag],na.rm=T),
  mean(dta22$WBPREJUDICE15[r.flag],na.rm=T)
)

#pdf("~/Dropbox/NAES/figures/prejudice-overtime-02182021.pdf")
plot(prejudice.scale.all,
     x=as.Date(c("2008-10-01","2012-10-01","2012-12-01","2016-02-01","2016-10-01","2016-12-01","2018-10-25","2020-01-31","2020-10-22")),
     type="b",ylim=c(-5,15),xaxt="n",pch=16,xlab="Date",ylab="Anti-Black Prejudice: Whites",cex.lab=1.3,
     xlim=as.Date(c("2008-01-01","2020-12-01")),yaxt="n",bty="n",lwd=2)
axis(side=2,at=seq(from=-5,to=15,by=5),labels=seq(from=-5,to=15,by=5))
axis(side=1,at=as.Date(c("2009-01-01","2010-01-01","2011-01-01","2012-01-01","2013-01-01","2014-01-01","2015-01-01","2016-01-01","2017-01-01","2018-01-01","2019-01-01","2020-01-01")),labels=c("2009","2010","2011","2012","2013","2014","2015","2016","2017","2018","2019","2020"), cex.axis=1.25)
par(new=T)
plot(prejudice.scale.d,x=as.Date(c("2008-10-01","2012-10-01","2012-12-01","2016-02-01","2016-10-01","2016-12-01","2018-10-25",
                                   "2020-01-31","2020-10-22")),
     type="b",ylim=c(-5,15),xaxt="n",pch=16,col="blue",yaxt="n",main="",xlab="",ylab="",xlim=as.Date(c("2008-01-01","2020-12-01")),bty="n",lwd=2)
par(new=T)
plot(prejudice.scale.r,x=as.Date(c("2008-10-01","2012-10-01","2012-12-01","2016-02-01","2016-10-01","2016-12-01","2018-10-25",
                                   "2020-01-31","2020-10-22")),
     type="b",ylim=c(-5,15),xaxt="n",pch=16,col="red",yaxt="n",main="",xlab="",ylab="",xlim=as.Date(c("2008-01-01","2020-12-01")),bty="n",lwd=2)
text(x=as.Date("2014-10-01"),y=prejudice.scale.r[3]+2,"Republicans",cex=1.2,col="red")
text(x=as.Date("2014-10-01"),y=prejudice.scale.d[3]-2,"Democrats",cex=1.2,col="blue")
#dev.off()

flag <- ! dta22$HELPBLK13 %in% c(NA) & ! dta22$HELPBLK15 %in% c(NA) & dta22$WHITE6==1
sum(flag)
d.flag <- ! dta22$HELPBLK13 %in% c(NA) & ! dta22$HELPBLK15 %in% c(NA) & dta22$PID4 %in% c(1,2,3)  & dta22$WHITE6==1
sum(d.flag)
r.flag <- ! dta22$HELPBLK13 %in% c(NA) & ! dta22$HELPBLK15 %in% c(NA) & dta22$PID4 %in% c(5,6,7)  & dta22$WHITE6==1
sum(r.flag)


helpblk.scale.all <- c(
  mean(dta22$HELPBLK6[flag],na.rm=T),
  mean(dta22$HELPBLK7[flag],na.rm=T),
  mean(dta22$HELPBLK13[flag],na.rm=T),
  mean(dta22$HELPBLK14[flag],na.rm=T),
  mean(dta22$HELPBLK15[flag],na.rm=T)
)

helpblk.scale.d <- c(
  mean(dta22$HELPBLK6[d.flag],na.rm=T),
  mean(dta22$HELPBLK7[d.flag],na.rm=T),
  mean(dta22$HELPBLK13[d.flag],na.rm=T),
  mean(dta22$HELPBLK14[d.flag],na.rm=T),
  mean(dta22$HELPBLK15[d.flag],na.rm=T)
)

helpblk.scale.r <- c(
  mean(dta22$HELPBLK6[r.flag],na.rm=T),
  mean(dta22$HELPBLK7[r.flag],na.rm=T),
  mean(dta22$HELPBLK13[r.flag],na.rm=T),
  mean(dta22$HELPBLK14[r.flag],na.rm=T),
  mean(dta22$HELPBLK15[r.flag],na.rm=T)
)

#pdf("~/Dropbox/NAES/figures/helpblk-overtime-02182021.pdf")
plot(helpblk.scale.all,
     x=as.Date(c("2012-10-01","2012-12-01","2018-10-25","2020-01-31","2020-10-22")),
     type="b",ylim=c(1,7),xaxt="n",pch=16,xlab="Date",ylab="Blacks Should Help Themselves",cex.lab=1.3,
     xlim=as.Date(c("2012-01-01","2020-12-01")),yaxt="n",bty="n",lwd=2)
axis(side=2,at=seq(from=1,to=7,by=1),labels=c("Gov't",seq(from=2,to=6,by=1),"Blacks"))

axis(side=1,at=as.Date(c("2012-01-01","2013-01-01","2014-01-01","2015-01-01","2016-01-01","2017-01-01","2018-01-01","2019-01-01","2020-01-01")),
     labels=c("2012","2013","2014","2015","2016","2017","2018","2019","2020"), cex.axis=1.25)
par(new=T)
plot(helpblk.scale.d,x=as.Date(c("2012-10-01","2012-12-01","2018-10-25",
                                 "2020-01-31","2020-10-22")),
     type="b",ylim=c(1,7),xaxt="n",pch=16,col="blue",yaxt="n",main="",xlab="",ylab="",xlim=as.Date(c("2012-01-01","2020-12-01")),bty="n",lwd=2)
par(new=T)
plot(helpblk.scale.r,x=as.Date(c("2012-10-01","2012-12-01","2018-10-25",
                                 "2020-01-31","2020-10-22")),
     type="b",ylim=c(1,7),xaxt="n",pch=16,col="red",yaxt="n",main="",xlab="",ylab="",xlim=as.Date(c("2012-01-01","2020-12-01")),bty="n",lwd=2)
text(x=as.Date("2014-10-01"),y=helpblk.scale.r[3]+.85,"Republicans",cex=1.2,col="red")
text(x=as.Date("2014-10-01"),y=helpblk.scale.d[3],"Democrats",cex=1.2,col="blue")
#dev.off()

#### regressions ----
### election fairness
dep.vars <- c(1,8,9,12,15)
out.ef <- list()
k <- 1
for(i in dep.vars){
  txt <- paste("dta23$DV <- dta23$ELECTFAIR",i,sep="")
  eval(parse(text=txt))
  out.ef[[k]] <- lm(DV ~ AGE6+EDYEARS6+INCOME6+FEMALE6+WHITE6+BLACK6+HISP6+
                      CATHOLIC0+PROTESTANT0+UNION0+PID6,data=dta23)
  k <- k+1
}
n.list <- length(out.ef)
texreg(list(out.ef[[1]],out.ef[[2]],out.ef[[3]],out.ef[[4]],out.ef[[5]]),stars=0.05,digits=3,
       custom.model.names=as.character(dep.vars))




idx <- which((! dta22$LEGIT6A %in% c(NA)) & (! dta22$LEGIT11A  %in% c(NA))  & 
               (! dta22$LEGIT15A  %in% c(NA))) 
length(idx)

dta23 <- dta22[idx,]
dep.vars <- c(6,8,9,11,15)
out.ld <- list()
k <- 1
for(i in dep.vars){
  txt <- paste("dta23$DV <- 6-dta23$LEGIT",i,"A",sep="")
  eval(parse(text=txt))
  out.ld[[k]] <- lm(DV ~ AGE6+EDYEARS6+INCOME6+FEMALE6+WHITE6+BLACK6+HISP6+
                      CATHOLIC0+PROTESTANT0+UNION0+PID6,data=dta23)
  k <- k+1
}
n.list <- length(out.ld)
texreg(list(out.ld[[1]],out.ld[[2]],out.ld[[3]],
            out.ld[[4]],out.ld[[5]]),stars=0.05,digits=3,
       custom.model.names = as.character(dep.vars))





idx <- which((! dta22$LEGIT6D %in% c(NA)) & (! dta22$LEGIT11D  %in% c(NA))  & 
               (! dta22$LEGIT15D  %in% c(NA))) 
length(idx)

dta23 <- dta22[idx,]

dep.vars <- c(6,8,9,11,15)
out.ld <- list()
k <- 1
for(i in dep.vars){
  txt <- paste("dta23$DV <- 6-dta23$LEGIT",i,"D",sep="")
  eval(parse(text=txt))
  out.ld[[k]] <- lm(DV ~ AGE6+EDYEARS6+INCOME6+FEMALE6+WHITE6+BLACK6+HISP6+
                      CATHOLIC0+PROTESTANT0+UNION0+PID6,data=dta23)
  k <- k+1
}
n.list <- length(out.ld)
texreg(list(out.ld[[1]],out.ld[[2]],out.ld[[3]],
            out.ld[[4]],out.ld[[5]]),stars=0.05,digits=3,
       custom.model.names = as.character(dep.vars))



dta22.sub <- dta22[! dta22$TRUMPFT10 %in% c(NA) & ! dta22$TRUMPFT15 %in% c(NA),]
dim(dta22.sub)
quantile(dta22.sub$TRUMPFT15,.9)

dep.vars <- c(10,11,13,14,15)
out.ft <- list()
k <- 1
for(i in dep.vars){
  txt <- paste("dta22.sub$DV <- dta22.sub$TRUMPFT",i,sep="")
  eval(parse(text=txt))
  out.ft[[k]] <- lm(DV ~ AGE6+EDYEARS6+INCOME6+FEMALE6+WHITE6+BLACK6+HISP6+
                      CATHOLIC0+PROTESTANT0+UNION0+PID6,data=dta22.sub)
  k <- k+1
}
library(texreg)
n.list <- length(out.ft)
texreg(list(out.ft[[1]],out.ft[[2]],out.ft[[3]],
            out.ft[[4]],out.ft[[5]]),stars=0.05,digits=3,
       custom.model.names = as.character(dep.vars))

### vote choice

dta23.sub <- dta22[! dta22$TRUMPVC11B %in% c(NA) & ! dta22$TRUMPVBIDEN15B %in% c(NA),]
dim(dta23.sub)

#dep.vars <- c(10,11,13,14,15)
vars <- c("ROMNEYVO6B","ROMNEYVO7B","TRUMPVC10B","TRUMPVC11B","TRUMPVC12B","TRUMPVBIDEN14B","TRUMPVBIDEN15B")

out.vt <- list()
k <- 1
for(i in 1:length(vars)){
  txt <- paste("dta23.sub$DV <- dta23.sub$",vars[i],sep="")
  eval(parse(text=txt))
  out.vt[[k]] <- lm(DV ~ AGE6+EDYEARS6+I(INCOME6/1000)+FEMALE6+WHITE6+BLACK6+HISP6+
                      CATHOLIC0+PROTESTANT0+UNION0+PID6,data=dta23.sub)
  k <- k+1
}
library(texreg)
texreg(list(out.vt[[1]],out.vt[[2]],out.vt[[3]],
            out.vt[[4]],out.vt[[5]],out.vt[[6]],out.vt[[7]]),stars=0.05,digits=3,
       custom.model.names = vars)

